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We put forward the use of total-variation- diminishing (or more generally, strong 



\^ stability preserving) implicit-explicit Runge-Kutta methods for the time integra- 

ls tion of the equations of motion associated with the semiconvection problem in the 

\^ simulation of stellar convection. The fully compressible Navier-Stokes equation, aug- 

mcntcd by continuity and total energy equations, and an equation of state describing 
^ the relation between the thermodynamic quantities, is semi-discretized in space by 

^ essentially non-oscillatory schemes and dissipative finite difference methods. It is 

•j-j subsequently integrated in time by Runge-Kutta methods which are constructed 

rS such as to preserve the total variation diminishing (or strong stability) property 

satisfied by the spatial discretization coupled with the forward Euler method. We 
analyse the stability, accuracy and dissipativity of the time integrators and demon- 
strate that the most successful methods yield a substantial gain in computational 
efficiency as compared to classical explicit Runge-Kutta methods. 
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Introduction 



Numerical hydrodynamical simulations are a common tool in astrophysical 
research. Just as some of their counterparts in the atmospheric sciences and 
in oceanography, astrophysical fluid flows are characterized by a vast range of 
timescales which are present in the solutions of the dynamical equations gov- 
erning the temporal evolution of such flows . Large relative changes of the 
solutions typically occur on the hydrodynamical timescale Tuuid = Ax/|u|. 
Here, u is the local flow velocity and Ax is the local spatial resolution of 
the simulation, which coincides with the grid size obtained from spatial dis- 
cretization of the governing partial differential equations. However, some of the 
physically important processes can also operate on much shorter timescales 
than Tiiuid- Examples include radiative transfer, sound waves, magnetohydro- 
dynamic processes, and chemical or nuclear reactions (see [23] for example, 
and references therein). 

In stellar astrophysics the two most important among those timescales are 
that of radiative energy exchange at the scale of a grid cell, Tj-ad, and the time 
Tsound a sound wave needs to cross such a cell. 

As long as sound waves are energetically or dynamically unimportant, a nu- 
merical simulation can be advanced with much larger time steps by using 
semi-implicit time integration methods, for example, by a fractional step ap- 
proach (see [H] for a general introduction). 

Similarly, if radiative transfer has the numerical characteristics of a stiff prob- 
lem, as is the case, for instance, for numerical simulations of the surface layers 
of A- type stars p9j, implicit time integration appears desirable as well. An- 
other important example where radiative diffusion can limit the time-step is 
the numerical simulation of double-diffusive processes in stellar interiors. 

Semiconvection is the most important special case of double- diffusive con- 
vection in astrophysics. Models of stellar structure and evolution predict set- 
tings where the heavier product of nuclear fusion provides stability to a zone 
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which otherwise would be unstable to convective overturning, because temper- 
ature sufficiently rapidly decreases against the direction of gravity. Such a zone 
would become convective if its composition were mixed. The question whether 
such a zone should be treated as if it were mixed or not is referred to as the 
semiconvection problem (see [26], [38], and Chap. 13.3 and 13-A in [1?], for 
example). A thorough physical analysis of the semiconvection problem based 
on numerical simulations in two spatial dimensions for a parameter set rele- 
vant to stellar astrophysics is given in |49] . Further discussions and reviews on 
this topic can be found, for instance, in [m5]l24|l43|H6] . For this problem, long 
total integration times are required even when the microphysics is idealised, 
whereas the time integration step is governed by Tj-ad and Tsound- 

For reasons outlined above, in this paper we discuss the advantages of implicit- 
explicit (IMEX) Runge-Kutta methods for simulations of stellar convection 
and diffusion in the parameter regime commonly associated with semiconvec- 
tion as discussed in [49] . These methods treat only part of the right-hand sides 
implicitly, where the resulting (generally nonlinear) equations can be solved by 
means of a generalized Poisson problem. It turns out that the total-variation- 
diminishing (TVD) property is essential for a numerical time integrator to be 
successful in simulations of the problems in our focus: to suppress spurious 
oscillations in the spatial discretization (which in this paper we realize for 
the hyperbolic terms by essentially non-oscillatory schemes and by dissipative 
centered finite difference schemes for the parabolic terms), this property has 
been demonstrated to be necessary for a stable integration in jl3j . 

The TVD property is more generally referred to as strong stability preserving 
(SSP) or monotonicity when norms other than the total variation norm or 
even sublinear functionals are considered. When the space discretization has 
the property that the functional of the discrete spatial profile is decreased 
in the course of numerical time propagation by the forward Euler method 
for a time-step AtpE, then an SSP method preserves this property under a 
step-size restriction of the form At < CAtpE with C > 0. Since the term total- 
variation-diminishing is more commonly used in the context of astrophysical 
simulations, where the total- variation-seminorm is the functional of interest, 
we mostly use these terms here synonymously. We expect the SSP (or TVD) 
IMEX methods to be useful also for other astrophysical problems where a high 
radiative (conductive) diffusivity of internal energy (temperature) restricts the 
time-step of hydrodynamical simulations, such as simulations of stellar surface 
convection with steep temperature gradients or at high resolution. 

The outline of the paper is as follows. First, we introduce SSP IMEX Runge- 
Kutta schemes and survey the related literature in Section [Tj 

Next, in Section|2]we specify the general set of equations to be solved in numer- 
ical simulations of semiconvection and related flows and describe the general 
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solution techniques implemented for this class of problems in the ANTARES 
code [HF] which we use for the numerical examples discussed further below. 
Subsequently, we discuss how this framework has to be modified when solving 
the dynamical equations with the IMEX approach. 

In Section [3] we analyse several SSP IMEX methods from the literature with 
respect to their radius of absolute monotonicity, stability and dissipativity. 
We show that the methods yield a significant advantage over classical explicit 
Runge-Kutta schemes with respect to both efficiency and accuracy and we also 
suggest a modification for one of the methods, which turns out to improve its 
efficiency. 

We then present numerical simulations of semiconvection in Section |4] to 
demonstrate the efficiency of the SSP IMEX methods as compared to the 
classical, explicit SSP Runge-Kutta time integrators and some non-SSP IMEX 
methods from the literature by giving numerical examples for a single layer in 
a physical scenario similar to that one studied in jlH]. We conclude this paper 
by a summary of the main properties of the IMEX methods, suggesting rea- 
sons for preferring particular methods and providing an outlook on interesting 
applications, which appear especially suited for this numerical approach. 



1 Implicit— Explicit Runge-Kutta Methods for Semiconvection 



To introduce our numerical methods in an abstract setting, we consider the 
ODE initial value problem 

yit)=F{y{t)) + G{y{t)), y{0) = yo, (1) 

where we assume that the vector fields F and G have different stiffness prop- 
erties. For this type of problems, partitioned Runge-Kutta schemes [15j, also 
called additive Runge-Kutta schemes, are popular. Methods of this kind use 
different Runge-Kutta formulae for the treatment of the two vector fields. We 



will see that the spatial semi-discretization of (15) below and the associated 



boundary conditions give rise to this kind of system. 

An s-stage partitioned Runge-Kutta method characterized by coefficient ma- 
trices A = (ttij) and A = {cLij) defines one step T/oid — ^ Z/now by 



yi = yo\A + At ^ aijF{yj) + At ^ aijGijjj), i = 1, . . . , s, (2) 

i=i i=i 

s s 

I/new = I/old + At 5] hF{y,) + At 5] hG{y,). (3) 
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If = for j > i, the method is referred to as an implicit-explicit (IMEX) 
method. 

These methods have first been investigated with respect to the SSP or TVD 
property in the context of hyperbohc systems with relaxation, where G = 
-G, £ ^ 1, in [36]. Ibidem, the common specification for strong stability 
preserving IMEX methods is introduced. An IMEX method is referred to 
as 'SSPfc(s, cr,p)' when it has the following properties: k is the order of the 
method in the stiff limit e 0, which is characterized by the coefficients 
for the explicit part. The latter must necessarily be SSP and is referred to 
as the asymptotically SSP scheme, s is the number of stages in the implicit 
scheme and a the number of stages in the explicit scheme, p is the global 
order of the resulting combined method. It is essential to observe that if the 
implicit scheme characterized by A = ( diagonally implicit Runge- 

Kutta (DIRK) method, then the explicit part is evaluated only once in each 
stage, providing the desired computational advantage 

The analysis in [36] is valid only for e ^ 1 |21]. However, several useful ex- 
amples of strong stability preserving IMEX Runge-Kutta methods are given, 
see Section IHl 

[21] develops a comprehensive theory of strong stability preserving additive 
Runge-Kutta schemes which extends the concepts for standard Runge-Kutta 
methods in a natural way: 

Let r, f be the step-size restrictions for monotonicity of the explicit Euler 
method for the vector fields F and G, respectively. We define the region of 
absolute monotonicity 

1Z{A,A) = {(r, f) G : (A, A) is absolutely monotonous on [— r, 0] x [— f, 0]}, 

(4) 

where the absolute monotonicity at a point (rg, tq) is characterized by al- 
gebraic relations for the matrices A, A. The boundary in the first quadrant, 
&R{A, A) n {(r, f) : r, f > 0}, is denoted as the curve of absolute monotonicity. 
The significance of the region TZ{A, A) is expressed in the following theorem 

m- 

Theorem 1.1 Let {A, A) be absolutely monotonous at (— r, — f) with step-size 
coefficients r, f. Then for h < min {rr,ff}, diminishing of the norm holds, 

WUiW < \\yo\d\\, i = l,...,S, Ibncwll < Iboldll- 

[22] gives order barriers for strong stability preserving additive Runge-Kutta 
methods similarly to [2H]- The order of an additive Runge-Kutta method 
{A, A) is bounded by the orders of A and A, respectively. This implies for 
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IMEX methods the order barrier p < 4 [28]. Moreover, [22] gives a simple 
algebraic criterion for a nontrivial region of absolute monotonicity in terms 
of incidence matrices of A, A. Some examples of strong stability preserving 
IMEX Runge-Kutta methods analysed in Section [s] can be seen in [221I28] . 

Some other relevant issues can be studied for IMEX methods. In this paper 
we focus on stability regions, error constants and dissipativity analysis. We 
close this section with a brief description of these concepts. 

The stability region of IMEX Runge-Kutta methods is defined in [T1I2] via the 
test equation of the form ([T|, where 

F{u) = i/3u, G{u) = au, a<0 < /3. (5) 

For this problem, 

Vnew = R{z)yoid, z = ttAt + i/3At, (6) 

and the stability region is the part of the complex plane where |-R(2;)| < 1. 
We will perform the corresponding analysis of the stability function for the 
methods considered in Section [HI 

To determine error constants of the methods we have computed the empirical 
convergence orders by solving the non-linear test problem 

y'{t) = (1 + sin(y(t))) + {y\t) - sm{y{t))), y{0) = 0, (7) 

with the known exact solution y{t) = tan(t). In this paper, the error constants 
are determined from the errors at t = 1.3. Their size is vital for the comparison 
of the accuracy of the methods of the same order and therefore the assessment 
of the work/precision relation. Of course, this single example only gives a 
rough indication of the size of the error constant and no rigorous estimate, 
but it seems sufficient for our purpose of comparing the methods in our focus 
with respect to accuracy, which we do in Sections |4| and |5] further below. The 
example was chosen such as to represent a nonlinear initial value problem with 
known solution whose profile is arbitrarily unsmooth as t — )■ |. 

Finally, we study the dissipativity of the time integrators in conjunction with 
suitable space discretizations. [45] gives a justification for considering only 
the diffusion term in this context since the advection term becomes negligible 
asymptotically. We will thus investigate the dissipativity of the implicit scheme 
specified by A. To this end, we apply the spatial discretizations in our 
focus to the heat equation 

Ut = bu^x, 

and associate for the spatial discretization Uj±k ^ e^'^'^. Thus, we compute the 
amplification factor g{9,^) = R{{LAxu)j), with n = h-r^^. This represents 
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the factor by which oscillations of frequency 9 are amplified in each time 
step. We pay particular attention to the case 9 = n corresponding to the 
mesh width. If = or |5f| = 1 also for a smaller < |^o| < ^r, then this value 
would represent the limit for a robust integration. However, such a pathological 
behaviour is only conceivable for methods of higher order. 

The spatial discretizations which were found to show a dissipative behaviour 
in |27j are the second-order three-point scheme (the upper index refers to the 
time step) 

{Axy 

and the fourth-order stencil 



Uxxi^-^ j : tn) ~ /A ^2 ' ("^^2;^ )j) 



(LaxU ), .= ^^^^ . (9) 



These are the methods actually implemented in ANTARES, where ^ is the 
default (see also [27]). 



2 Solving the Hydrodynamical Equations with ANTARES 



In our model, the fundamental equation of motion is the fully compressible 
Navier-Stokes equation which describes momentum conservation: 

(pu)' + V ■ (pu ® u + pJ) = pg + V • a. (10) 

The state variables in the model equations generally depend on the spatial 
variables (x, y, z) and time t. In the simulations presented in Section |4] we 
solve problems in two spatial variables, whence the variable z will be dropped 
in the rest of the paper. The (explicit) dependencies are stated in Table [l| For 
simplicity, we omit the dependencies in the problem specification ([Io|)-([l4)). 
The model is completed by the continuity equation 

P+V-(pu) = 0, (11) 

which ensures conservation of mass, and the total energy equation 

e' + V ■ (u(e + p)) = p(g • u) + V • (u • (t) + g^ad, (12) 

which describes conservation of the latter. In the case of a two-component 
fluid, the system is augmented by the concentration equation of the second 
species, 

(cp)' + V- (cpu) = V- (p/tcVc). (13) 
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The variables and parameters which appear in the model formulation are 
collected in Table [U 



p = 


p{x,y,z,t) 


gas density 


c = 


c{x,y,z,t) 


concentration of second species 


u = 


u{x,y,z,t) = {u.v.wY 


flow velocity 


pu 




momentum density 




u 


Kronecker product 


p = 


P{T, p) 


gas pressure 


g = 


(5,0,0)^ 


gravitational acceleration 


a = 


a{x,y,z,t) 


viscous stress tensor for zero bulk viscosity 


V 




dynamic viscosity (appears in the definition of a) 


e = 


e(x, y, z, t) = eint + ekin 


total energy density 


T = 


: T{x,y,z,t) 


temperature 


Qrad Qrad(''^) ?/; ^i ^) 


radiative source term 


Cp = 


= Cp{T,p,c) 


specific heat at constant pressure 


Xi^ -- 


= Xu{T,p,c) 


(specific) opacity at frequency v 


K = 


--K{T,p,c) 


radiative (or thermal) conductivity 


Kt 


= K/ic^p) 


radiative (or thermal) diffusivity 


Kc - 


= 'ic{T,p,c) 


diffusion coefficient for species c 


Iu = 


= 4(r), r = r{x,y,z) 


specific intensity along the ray of direction r 


Su - 


= Su{x,y,z) 


source function 



Table 1 

Variables and parameters in the equations ((lO|)-( 14 ) 



In general, the radiative source term Qrad is determined as the stationary limit 
of the radiative transfer equation 



r • Vh = PXv{Sy - lu), 



(14) 



which is solved for all ray directions r and for all frequencies u, resulting in the 
specific intensity /^(r), for details see |17j. Si, here denotes the source function. 



The equations of hydrodynamics ( 10 ), ( 11 ) and ( 12 ) are closed by the equation 



of state which describes the relation between the thermodynamic quantities. 
For the particular choice, see [35] . 

For the initial condition, a slightly perturbed static model stellar atmosphere 
or stellar envelope is used which is equipped with a small seed velocity field 
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or density perturbation to start dynamics away from equilibrium. 



In the framework discussed below, boundary conditions are based on the as- 
sumption that all quantities are periodic in both horizontal directions. More- 
over, for the hydrodynamical equations, "closed" (Dirichlet) boundary condi- 
tions at the upper and lower boundary of the computational domain are used. 
A recent development is to replace these by "open" (Robin) boundary condi- 
tions. These allow inflow and outflow of fluid along the vertical direction which 
is defined by the direction of g. For the radiative transfer equation (14), incom- 
ing radiation at the boundary of the computational domain must be specified. 
Since double-diffusive convection in stars takes place in regions which are op- 
tically thick, the quantity Qrad can accurately be obtained by means of the 
diffusion approximation for radiative transfer, Qrs,A = V ■ Frad = V ■ {KVT). 
In this case, further knowledge about the intensity ly is not necessary. 



The ANTARES code [35] solves this system of equations numerically in either 
one, two, or three spatial dimensions on a rectangular grid (spherical coordi- 
nates with a logarithmically rectangular grid are also possible, i.e., the grid 
may be locally rectangular with logarithmic grading in the radial component). 



For the spatial discretization, ANTARES allows the definition of several grids 
which can be nested inside each other to improve resolution in regions of 
interest. At the moment, ANTARES provides up to three levels of nested grids. 
For the hyperbolic terms, discretizations of ENO {essentially non- oscillatory 
PU] ) type are implemented. These comprise classical ENO methods, WENO 
[weighted essentially non- oscillatory) methods |10] (optionally in conjunction 
with Marquina flux splitting [8]) and CNO {convex non- oscillatory) schemes 
|31j . Each of the methods uses adaptive stencils which are chosen such as to 
avoid spurious oscillations in the computed solution. The spatial derivatives 
are calculated for each direction separately. 

The parabolic terms are discretized by dissipative finite difference schemes [27] 
of fourth order. The radiative heating rate is determined by the short char- 
acteristics method, or by means of the diffusion approximation for radiative 
transfer, where appropriate. For the time integration, total variation dimin- 
ishing Runge-Kutta methods [3911^ are employed. 

ANTARES implements two different parallelization concepts. For architec- 
tures with distributed memory, domain decomposition is used and realized by 
an MPI implementation. In this approach each grid is split along the horizon- 
tal direction(s) and optionally, also along vertical ones, into subdomains. The 
memory required to store the computational variables for each subdomain is 
provided by the resources available to the CPU core performing the computa- 
tions necessary for that subdomain. In this way, each CPU core is mapped to 
a specific geometrical volume. However, since some supercomputers offer only 
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a limited amount of memory per CPU core and because the domain decompo- 
sition approach is not very efficient on small grids, ANTARES offers a second 
type of parallelization which can be used along with or independently of the 
former. It is based on a shared memory concept for each subdomain and is 
implemented through OpenMP directives. Therefore, the most time consum- 
ing operations which can also be performed independently of each other are 
identified and computed in parallel. This approach scales only to a moderate 
number of CPU cores, but allows improvement of the scaling and the compu- 
tational speed of the domain decomposition based parallelization for a larger 
number of problems and for a greater variety of computer architectures. 



In the following, the dynamical evolution of the fluid is described by the mul- 
tispecies Navier-Stokes equations presented above. Additionally, dimension- 
less quantities such as the Prandtl number Pr = c^rijK^ the Lewis number 
Le = CppHc/K, the Rayleigh number Ra and the stability parameter Rp are 
defined to determine the diffusivities kt-, k,c and the viscosity rj. The former 
quantities arise in the definition of the starting model but do not appear in 
the evolution equations ( 15 ) below. Since we solve the dynamical equations for 
a compressible flow, we specify the vertical extent of the simulation domain 
in multiples of the pressure scale height Hp = P/{pg). For the simulations 
presented in Section [4] the domain always covers IHp. In the physical model, 
intermolecular forces are neglected, so the fluid is assumed to be an ideal gas. 
The radiative source term Qrad is modelled using the diffusion approximation 
with a heat conductivity K which is constant in time and otherwise only a 
function of the vertical coordinate [33fM] . We use a variant of this setting 
here, where not only Pr, Le, and Cp, but also K, Kc, and rj/p take constant 
values [l9]. This setup simplifies studies of the basic physics while it is still 
useful for extrapolations to astrophysically relevant cases (cf. \49\ ) . 



For the model problem, the multispecies Navier-Stokes equations can be recast 
as 



d 
dt 



pc 
pu 

V ^ J 



-V- 



pu 
pcu 

pu (g) U + P — (T 

y eu + Pu — w ■ a I 



( \ 


99 
\P9^ ) 



V 





pK(S7c 


K\/T 



Giy{t)) 



(15) 



In the context of the problem (|15|), the i^^ implicit stage of an IMEX Runge 



Kutta method is typically of the form 
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yi=y* + Atai^iG{yi), 

where y* is known from previous stages. This is a consequence of aij 
j > i. 

This translates to 



(16) 
for 



Pi = P , 

(pc). = {pc)* + Atdi^iV ■ {piKcVci 

{p\i)i = (pu)*, 

e, = e* + Atai,,V-(irVT,). 



(17) 
(18) 
(19) 
(20) 



Rearranging (18) leads to 



At a 



— Ci-V- {p*KcVCi 



p'^e 

At cij 



(21) 



Obviously, this is a general elliptic equation for q of the form 



g{x, y)ip{x, y)-V ■ {h{x, y)V(fi{x, y)) = f{x, y). 



(22) 



Due to model assumptions, (20) can also be transformed to resemble a general 



elliptic equation. We start by recalling 



6 ^int ~l~ *^kin 

1 , , 

= eint + oP|u| 



(23) 
(24) 



Bearing in mind equation (19), equation (20) reads 



e*„t + AtS,,,V-(KVr,). 



(25) 



The equation of state for an ideal gaQ relates the temperature T and the 
internal energy Cint via 

_ 3 Tp_Rgas 



2 m 



(26) 



^ Note that, as is common in astrophysics, the gas constant Rgas is taken to be 
relative to the atomic mass unit such that the dimensionless mean molecular weight 
can be used in the equation of state instead of the molar mass. 
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if we assume the ratio of the specific heats at constant pressure and volume 
to equal 5/3. Here, m denotes the mean molecular weight of the compound. 



So we have at stage i 



eint. = !^^^^ (27) 
2 rrij 



and therefore, we arrive at 



-RgasP 



Since rrii is evaluated using the mass fraction Cj, it is necessary to solve equa- 



tion (21) first 



Thus, the solution of an implicit stage translates to the solution of the gener- 
alized Poisson equations for the mass fraction c and the temperature T. In the 



ANTARES framework, finite elements are used for the discretization of (22). 
The resulting linear system is solved by the conjugate gradient method. For 
parallel computations, the Schur complement algorithm is applied. A detailed 
description is given in |14j . 

The above procedure applies without modification to the case of the fully 
compressible Navier-Stokes equation. However, for low Mach number flows 
a splitting approach is preferable where the terms containing pressure are 
treated separately in a post-processing step, i. e. after evaluation of all other 
terms for the computation of the velocity fields. The latter are obtained for the 
current step from an additional generalized Poisson equation for the pressure. 
This procedure is described in detail in [TT]. Consequently, the explicit stage 
is evaluated here using a fractional step method [30] as implemented in [T7] . 



3 Strong Stability Preserving IMEX Schemes from the Literature 



In this section we study different SSP IMEX methods from the literature 
focusing on the topics described in Section [TJ The results are summarized in 



Figure 16 and Table 10 



All the strong stability preserving IMEX schemes listed in the following subsec- 
tions have DIRK (diagonally implicit Runge-Kutta) methods as the implicit 
scheme. This structure ensures that the stages can be solved successively, and 
the explicit part only has to be evaluated once in each stage. 
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An SSP2(2,2,2) Method 



gives an IMEX SSP2(2,2,2) method with nontrivial region of absolute 
monotonicity (7 = 1 ~ ;^)- 



A 





1 



7 7 
1 — 7 1 — 27 7 



(29) 



1 1 

2 2 



A 



1 1 

2 2 



The coefficients imply 7^(A) = 1, n{A) = 1 + ^2, and 

n{A, i) = {(r, f) : < r < 1, < f < v^(l - r)}, 



see 



The stability region is entirely located in the left half plane, tangent to the 
imaginary axis and unbounded as 3fJ(z) — >■ —00. Hence, the schemes are A{a)- 
stable with a = |, but not A-stable [16]. Moreover, limsR(2)_>._oo -R(-2) = 0. A 
plot of the stability region is given in Figure [T] (left). 





Fig. 1. Stability regions of IMEX method (29) (left) and A (right) 



The stability function of the implicit scheme A is 

(1 + v^) (1 + ^ + V2' 



Ra{^) = 2 



(2-Z + V2 



(30) 



A plot of the related stability region is shown in Figure [T] (right). The scheme 
A appears to be A-stable and satisfies lims({(2)^_oo R^i^) = 0, implying L- 
st ability. 

The dissipativity analysis for the implicit scheme defined by A yields the am- 
plification factors for the standard three-point space discretization ([8| and the 
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Table 2 



e 







1 


7r 


^ (l+V2)(l-2/j+^^+v^) 


4 


(-2-2/i+/x V2-v^)^ 


7r 
2 


(1+v^)(1-2m+v^) 
(2+2^+^2) 


vr 


^ (1+v^)(1-4m+v^) 
(2+4 M+V2)^ 



laDie / 

Values oi g{9, //) for some 9, implicit scheme in (29 ), three point space discretization 
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1 


TT 

4 


(1+^2) (6-15 M+8v^M+6V2) 


(-12-15 At+8 V2)^ 


TT 

2 


g (1+^2) (3-7^+3^2) 
(6+7^+3^2)^ 


vr 


(1+^2) (3-16 At+3V2) 
(6+16^4+3^2) 



Table 3 

Values of g{9, fj.) for some 9, implicit scheme in (29 ), fourth order space discretization 



fourth order stencil (|9]), respectively. The amplification factors are evaluated 
at the points 9 G {0, |, |,vr} in Tables |2] and [sj respectively. 

The first positive zero of g{n,n) is ~ 0.6035 for the three-point scheme (|8|, 
where the function changes its sign, and |(7(vr, never exceeds 1. The first 
positive zero of (^(vr, /i) for the fourth order space discretization (|9]) is ~ 0.4526, 
where the function changes its sign. The modulus never exceeds 1. 



Modification 0/7 



We may conceive of optimizing the method (29) by adapting the value of the 
parameter 7 in the definition of A according to the resulting stability, accuracy, 
and dissipativity properties. The region of absolute monotonicity depends on 
7 as follows 



14 



7^(i) 



47-1 



1 

1-37' 
1-27 

272-47+1 Y (272-47+1)^ 
1-27 



0<7< i 

1/4 < 7 < 1 - ^, 



7 = 1 



V2' 



272-47+1 



+ 



47-1 



(272-47+1)^ 



h<i<h 



V2 



n{A,A) 



{(r,f) : < r < 1 ,0 < f < If^} , < 7 < |, 



{ (r, f) : < r < 1^ , < f < ^ } , | < 7 < I- 
A plot of the function TZ{A) in dependence of 7 is given in Figure M 



Radius of absolute monotonicity 



3.5 
3 

« 2.5 
2 
1.5 

1 



0.1 0.2 0.3 0.4 0.5 

gamma 



(31) 



(32) 



Fig. 2. Radius of absolute monotonicity TZ{A) as a function of 7 for (29) 



The regions of absolute monotonicity 71{A, A) for the values 7 G {0.1, 0.2, 0.3} 
are plotted in Figure |3) 




0.2 0.4 0.6 0. 




Region of ab^oliuc 



0.2 0.4 0.6 0.8 1 




0.2 0.4 0.6 0.8 1 



Fig. 3. Regions n{A,A) for 7 E {0.1,0.2,0.3} for (29). 



The stability regions for the IMEX schemes for the different values of 7 cover 
bounded subdomains of the left half plane for 7 < 0.25, while for 7 > 0.25, 
the stability regions cover unbounded domains in the left half plane. In fact, 
the left boundaries zieh satisfy 



^^left 



47-1 ' 



7 < 0.25, 



-00, 7 > 0.25. 
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However, even in the cases with unbounded stabihty regions, in general 
hm: 



3 R{z) 7^ 0. The boundaries of the stabihty regions are plotted in 
Figure |4] (left), where values equal to represent unbounded stability regions. 
The implicit schemes A show the same stability behaviour concerning both 
the boundaries of the stability regions and the limits for ^{z) — ?► — oo. 





-10 
-20 
-30 
^0 
-50 
-60 



Left boundaiy of stability region 



Error constant 



0.2 0.3 0.4 0.5 



0.3 
gamma 



Fig. 4. Left boundary of the stability region (left) and error constant computed for 
^ (right) as a function of 7 for (29). 



It was demonstrated with a Matlab implementation that the convergence 
order two is retained also for the modified values of 7. The error constant 
depends on 7, however. In Figure |4] (right) we plot the error constant as a 
function of 7, where the error is determined at t = 1.3 for the test problem 
([7]). We note that for small 7, the error constant decreases as 7 grows, while 
for 7 > 0.1833 the constant grows monotonically. This behaviour does not 
appear to be related to the results we obtain for the dissipativity analysis, 
where for 7 = 0.25 the behaviour changes. 



Roots and maxima 



16 
14 
12 
10 
1 8 
6 
4 
2 



roots 
maxima 



0.3 



0.1 0,1s 0.2 



Fig. 5. Dissipativity analysis for (29) for different values of 7, three point space 



discretization ([s]) (left) and fourth order discretization ^ (right). 

To assess the dissipativity properties of the modified scheme we vary 7 and 
compute the first positive zero of g{TT,fi) and the point /x where the modulus 
of this function exceeds 1. 



For 7 > 0.25, the amplification factor for the three-point space discretization 
(js) has a zero at /i = ^g^2liQ^'!^ , but there is no real root for 7 < 0.25. Also 
for 7 < 0.25, \g\ exceeds 1 at = ^3^, and this behaviour does not occur for 
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7 > 0.25. For the discretization (|9]), the behaviour is the same, but the scale 
with respect to /i is multiphed by 0.75 (see also Section 3.2 in [2Z])- 

To illustrate this analysis, in Figure [5] we vary 7 and compute for 30 values of 
7 spaced equidistantly in the interval [0.15, 1.1 — l/-\/2] the first positive zero 
of (7(7?, fi) and the point fi where the modulus of this function exceeds 1. This 
analysis is given for the three point space discretization ([s]) in Figure [s] We 
conclude that it may be of interest to choose 7 < 0.25 to ensure < (/(vr, fi) < 
1. This can be achieved by a value of 7 just slightly smaller than 0.25. 

Taking into account the results displayed in Fig. |4] we suggest to optimize 7 
by taking it as large as necessary to avoid any linear stability restrictions due 
to the term G{y{t)) in ([T| and as small as possible to minimize the stability 
constant C. Note that some special values for 7 are 0, (l-l/v^)/2 ^ 0.2113, 
and 1/4, in addition to the originally proposed value of 1 — l/\/2 ^ 0.2929. 
They are readily identified to yield the classical explicit SSPRK(2,2) methocp] 
of order two by Shu & Osher [51], the optimal implicit third order SSP method 
with two stages [ID], and the optimal implicit second order SSP method with 
two stages [lO] . In our numerical tests reported in Section |4] below, we found 
that the choice 7 = 0.24 yielded the most efficient time integrator. 



An SSP2(3,3,2) Method 



pi] gives an IMEX SSP2(3,3,2) method with nontrivial region of absolute 
monotonicity: 









1 

5 


1 
5 





1 

2 


iOO 


3 
10 


1 

10 


10 


1 


1 i 

2 2^ 


1 


1 

3 


1 1 

3 3 


A 


111 
3 3 3 


A 


1 
3 


1 1 
3 3 



This is a modification of a scheme from [36], where the latter turned out to 
have a trivial region of absolute monotonicity. It holds that 'R-{A) = 2 and 
R{A) = 1(^70-4), and 

n{A, i) = {(r, f) : < r < 1, < f < 0(r)}, 

where 

0(r) = |^(-28 + 9r) + ^ Vl264 - 984r + 201r2| . 

^ We will henceforth use the common specification 'SSPRK(s,p)' introduced in [28] 
for an s-stage order p explicit strong stability preserving Runge-Kutta method. 
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We note that the latter is a correction with respect to [21], since we have 
found r to be necessarily bound by 1 in TZ{A, A). A plot of 7^(A, A) is given 
in Figure |6| 



Region of absolute monotonicity 




Fig. 6. Region of absolute monotonicity for method (33) 



A plot of the stability region is given in Figure [7] (left). We observe that the 
stability region is tangent to the imaginary axis and appears to be unbounded 
as ^{z) — 7- — oo. Moreover, limsj{(2)^_oo R{z) = 0. The stability function of the 




Fig. 7. Stability region of IMEX method ([SSj) (left) and A (right), 
implicit scheme is given by 



R^{z) 



-5 + zY(-3 + z) 



(34) 



A plot is shown in Figure [T] (right). Note that the stability region is not 
connected. The method is A-stable, however. Moreover, the scheme A satisfies 



limi 



«(2)- 



0. 



The dissipativity analysis for the implicit scheme defined by A yields the 
amphfication factors for the standard three-point space discretization ^ and 
the fourth order stencil (|9]). These amplification factors are evaluated at the 
points 6* G {0, |, |, vr} in Tables 111 and Isl respectively. 



The first positive zero of g{TT,fi) is ~ 0.6064 for the three-point space dis- 
cretization ([s]), where the function changes its sign, and \g{Tr,fi)\ never ex- 
ceeds 1. The first positive zero of g{n,^) is ~ 0.4551 for the fourth order 
space discretization ([9]), where the function changes its sign, and \g{Tr,fx)\ 
never exceeds 1. 



18 



Table 4 
Values of g{ 



e 


9(0, f^) 





1 


TV 

4 


75+20 M V2-40 /i-27 ^^+18 ^t^v^ 


(-5+^^2-2 Aj)^(-3+At V2-2n) 


TV 

2 


-75+18^^+40^4 
(5+2 Aj)^ (3+2/^) 


vr 


-75+80 M+72 
(5+4aj)2(3+4/x) 



//) for some 0, implicit scheme in (33), three point space discretization 



9 







1 


TV 


g 1800-1200 /i+640/iV2-1059At2+720Ai2^ 


4 


(-30-15 M+8 At v^)^ (-18-15 /i+8 At v^) 


TT 

2 


n /o -450+280 At+147At^ 
/ (7At+15)^(7At+9) 


TT 


n 320At+384At2-225 
(15+16 At)^{9+16 At) 



Table 5 

Values of g{9, fi) for some 9, implicit scheme in (33 ), fourth order space discretization 
@- 

An SSP3(3,3,3) Method 



gives the following SSP3 (3,3,3) method with nontrivial region of absolute 
monotonicity: 



A 





1 



1 1 

4 4 







112 
6 6 3 
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14 


1 







15 


15 


1 


7 


1 


1 


2 


30 


5 


15 


A 


1 


1 


2 


6 


6 


3 



It holds TZ{A) = 1 and 7^(i) = ^(13 - 2^7), and 



where 



n{A, A) = {(r, f) : < r < 1, < f < 0(r)} 
15 



(/)(r) 



302 



28 - 25r - Vl80 - 192r + 21r2) . 



(35) 



A plot of 'R{A, A) is given in Figure [8 

The stability region occupies a bounded domain in the negative half-plane, 
and slightly overlaps the imaginary axis, see Figure M (left). Note that 



19 






Fig. 9. Stability regions of IMEX method (35) (left) and implicit method A (right) 



limsR(2)_!._oo = oo. The stability function of the imphcit scheme is given 

by 



450 + 390 2 + 167^2 _^ 47^3 
2 (-15 + 2)^ 



(36) 



A plot is shown in Figure [9] (right). The point where the stability region 
intersects the negative real half-line is located at x ~ —3.248 for the IMEX 
scheme. The same value is computed for the stability region of the implicit 
scheme. However, the stability region of A extends further along the imaginary 
axis. Finally, lim^(^^)^_^ l-R^l^)! = oo. 



The dissipativity analysis for the implicit scheme defined by A yields the 
amplification factors for the standard three-point space discretization ^ and 
for the fourth order stencil These amplification factors are evaluated at 
the points 6' G {0, |, |, vr} in Tables [g] and [7| respectively. 



For the three-point space discretization ([8]), the first positive zero of g{'n',n) 
is ~ 0.4650, where the function changes its sign, and g{7i,fi) = —1 for fi ^ 
0.8122. The first positive zero of g{7r, fi) is ~ 0.3488 for the fourth-order space 
discretization ([9]), where the function changes its sign, and (?(7r,/i) = —1 at 
/i = 0.6093. 
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vr 



225+195 M V2~390 fi+501 ^^-334 /j^ \/2+329 /j^ v^-470 /^^ 



225+390 At-334 ^j2 + 188 
(15+2^)^ 



1504^^+780^-1336 At^-225 
\ (15+4^)^ 



Table 6 

Values ol g{0, fj,) for some 9, implicit scheme in (35), three point space discretization. 



9(0, f^) 



1 



-^1^2 97200-210600 ^+112320 ^ \/2+353706 ^^-240480 /j^\/2-429345 At^+301928 ^l^V2 

(-90-15^+8^^2)^ 
- 12150+24570 /j-24549 At^+16121 /j'^ 



-1/6 



(45+7 At)" 



vr 



-1/3 



-6075+28080 /j-64128 /i^+96256 /t^ 
(45+16 m)^ 



Table 7 

Values oi g{6,fj,) for some 6*, implicit scheme in (35), fourth order space discretiza- 
tion. 



4 Numerical Experiments 



In this section, we present the results of the simulations performed with the 
time integrators discussed in this paper and compare their performance to the 
classical explicit 2-stages second order and 3-stages third order methods of 
Shu & Osher |?T] , and the 3-stages second order method of [2H] ■ They are the 
optimum SSP RK methods for their given number of stages and order and 
we refer to them by their common technical specifications 'SSPRK(2,2)' and 
'SSPRK(3,3)', as well as 'SSPRK(3,2)', respectively. The coefficients of the 
SSPRK(2,2) and SSPRK(3,3) methods were derived for a different purpose in 
[18j and in [U], where the third order method was proposed as an embedding 
formula for the second order method (see also [3]). Shu & Osher [H] derived 
a first framework for deducing higher order Runge-Kutta methods with the 
total variation diminishing property and first identified the optimal explicit 
second and third order methods with two and three stages, respectively. For 
these SSPRK(2,2) and SSPRK(3,3) schemes an analysis of their stability, dis- 
sipativity and accuracy properties can be found in [27]. For the SSPRK(3,2) 
method, a similar study is given in the Appendix section further below. Note 
that SSPRK(2,2), SSPRK(3,2), and SSPRK(3,3) are the exphcit schemes in 



the IMEX methods (29), (\33L and (35), respectively. 



Furthermore, we also consider some non-SSP methods from the literature. 
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both explicit and IMEX, for the sake of comparison, and also an asymptotically 
stable SSP IMEX scheme from 



Implementation Issues 

To define the test problem, according to [33j we specify a hydrostatic con- 
figuration which is unstable against convection. The simulation of a single 
semiconvective layer requires the mean molecular weight to be linearly (and 
stably) stratified. As time evolves, we expect convection to set in and mix 
the zone completely, although its development is inhibited by the stable mean 
molecular weight gradient. A critical quantity in this process is hence the 
buoyancy timescale and we return to this topic further below. 

The simulations shown here have been performed on the Vienna Scientific 
Cluster, using 64 CPU cores in parallel. The spatial resolution is 400x400 grid 
points. Simulation time is measured in units of sound crossing times (sort). 
One scrt is defined as the time taking an acoustic wave to propagate from the 
bottom to the top of the simulated box. In our simulations, 1 scrt = 5215.5 s 
and the simulation time is 200 scrt. 

Restrictions on the time-step At are imposed by heat diffusion tt, diffusion 
of the second species Tc, the viscosity Tvisc and the velocity of the fluid Tfiuid, 

At = min{rc, r^, r^isc, ra^id}, (37) 

where 

Te = — min{(Ax)^ {Ayf}, tt = — min{(Ax)^ {Ay)'}, 

TVisc = min{(Ax)^ (Ay)'}, Tfiuid = vTTx min{Aa;, Ay}, 

1] max(|u|) 

with (not necessarily equal) Courant numbers (CFL numbers) Cc, Ct, Cvisc, 
Cfluid- This assumes that the time-step limitation due to sound waves has been 
removed by a fractional step approach as mentioned at the end of Section |2] 
Otherwise, Tfluid additionally depends on the sound speed. We note that the 



source term in F{y{t)) in (15), which represents buoyancy forces acting on 
the flow, can be neglected in the limit where max {Ax, Ay} — )■ 0, since its 
contributions are of lower order (see I45j for a discussion of the treatment of 



lower order terms in stability analyses). 



Due to ([15]), IMEX methods treat the terms V ■ (pK^Vc) and V ■ (KVT) 
implicitly, so the restrictions and do not have to hold. Since at least the 
flrst part of simulations of semiconvection is usually dominated by diffusion 
processes and the Prandtl and Lewis numbers satisfy Pr < 1, Le < 1 in a 
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stellar context, the simulations are initially restricted by tt- Hence, IMEX 
methods can provide the desired computational advantage. 

To enhance the stability and the efficiency of our methods we have imple- 
mented a heuristic to adaptively select the time steps. The most effective 
criterion for regulating the time steps turned out to be to monitor two-point 
instabilities appearing in the conservative variables {p, pc, pu, Et). Due to the 
gravitational force operating vertically, such instabilities are prone to appear 
in the horizontal direction. 

To detect the occurrence of such oscillations, for each variable we use the 
difference between two grid cells to determine the sign of the corresponding 
gradient. In case of the density p this reads 

dl = Pi,j-1 — Pi,j-2, 
^■2 = Pi,j ~ Pi,j-li 

ds = Pi,j+i - Pi,j, 

di = Pi,j+2 ~ Pi,j+l) 



where 1 < i < Ux and 1 < j < Uy, assuming the grid consists of x Uy points. 
If the sign pattern of (di, d2, d^, ^4) corresponds to (+, — , +, ±), (— , +, — , ±), 
(±,+,— ,+,) or (±,— ,+,—), we have located a two-point instability. Since 
the time-step is initially chosen to be that one required for a fully explicit 



time integration method as in (37), such patterns are smoothed out rapidly, if 
present in the initial condition. Consequently, their later occurrence is a good 
indicator for an instability developing because of too large a time-step taken 
during the time integration. 

The time-step control permits the occurrence of n^^ ■ 0.1 two-point instabilities 
for fixed i. If this limit is exceeded, the time-step is repeated using a step-size 
decreased by a factor |. 

To permit the system to readjust after reducing the time-step no modifications 
of At are made for the next 15 time-steps regardless of the number of two- 
point instabilities. If the number of oscillations still exceeds the given limit 
after those 15 time steps, the time-step is again reduced. If no or very few 
two-point oscillations are encountered for over 50 successive time-steps, the 
time-step is augmented by a factor of |. 

Alternatively to this heuristic control, we could also monitor the rate of change 
in the solution to adjust the time-steps. However, the proper rate required to 
prevent the development of two-point instabilities for the present application 
turned out to be too pessimistic to achieve CFL numbers as high as discussed 
below. Furthermore, a simple control of the residual did not yield a satisfactory 
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Method 




^^mean 


CFLmax 


CFLmean 


CFLstart 


Singlelayer Pr = 0.1, Le = 0.1, i?p = 1.1, Ra* = 160000 


SSPRK(2,2) 


3.71 s 


3.71 s 


0.2 


0.2 


0.2 


SSPRK(3,2) 


9.31 s 


9.31 s 


0.5 


0.5 


0.5 


IMEX SSP2(2.2.2) 


22.21 s 


11. 5G s 


1.20 


0.62 


0.2 


IMEX SSP2(2,2,2), 7=0.24 


36.35 s 


19.44 s 


1.96 


1.05 


0.2 


IMEX SSP2(2,2,2), 7=0.24 


36.35 s 


19.43 s 


1.96 


1.05 


0.3 


IMEX SSP2(2,2,2), 7=0.24 


37.02 s 


19.45 s 


2.00 


1.05 


0.4 


IMEX SSP2(2,2,2), 7=0.24 


35.53 s 


19.47 s 


1.91 


1.05 


0.5 


IMEX SSP2(3,3,2) 


74.52 s 


57.37 s 


4.02 


3.09 


0.4 


IMEX SSP2(3,3,2) 


93.15 s 


57.16 s 


5.02 


3.08 


0.5 


SSPRK(3,3) 


3.71 s 


3.71 s 


0.2 


0.2 


0.2 


IMEX SSP3(3,3,3) 


15.14 s 


10.14 s 


0.82 


0.55 


0.2 


Singlelayer Pr = 0.5, Le = 0.1, Rp = 1.1, Ra* = 160000 


SSPRK(2,2) 


3.72 s 


3.72 s 


0.2 


0.2 


0.2 


SSPRK(3,2) 


9.31 s 


9.31 s 


0.5 


0.5 


0.5 


IMEX SSP2(2,2,2) 


23.14 s 


13.84 s 


1.24 


0.74 


0.2 


IMEX SSP2(2,2,2), 7=0.24 


23.13 s 


15.72 s 


1.24 


0.85 


0.2 


IMEX SSP2(2,2,2), 7=0.24 


22.76 s 


15.79 s 


1.22 


0.85 


0.3 


IMEX SSP2(2,2,2), 7=0.24 


20.32 s 


15.81 s 


1.09 


0.85 


0.4 


IMEX SSP2(2,2,2), 7=0.24 


22.81 s 


15.85 s 


1.23 


0.84 


0.5 


IMEX SSP2(3,3,2) 


40.65 s 


33.54 s 


2.19 


1.80 


0.4 


IMEX SSP2(3,3,2) 


40.65 s 


33.87 s 


2.19 


1.82 


0.5 


SSPRK(3,3) 


3.72 s 


3.72 s 


0.2 


0.2 


0.2 


IMEX SSP3(3,3,3) 


15.05 s 


9.70 s 


0.81 


0.52 


0.2 



Table 8 

Comparisons of time steps and CFL-numbers over the first 80 scrt for the case where 
Pr = 0.1 and over the entire 200 scrt for the case where Pr = 0.5 (see also text). 



behaviour. This leaves the heuristic time-step control as the most effective 
method for the IMEX based time integration of our simulations. 



24 



Numerical Results with SSP Schemes 



Tables |8] and |9] sum up the performance achieved with the presented Runge- 
Kutta schemes. The evolution of the time steps is illustrated graphically in 



Figures [TO] and [TT] for the two simulation scenarios we have focused on. Since 
the capability of the method is best judged in that part of the simulation 
where the fluid velocity is too small to severely limit the time-step At, the 
CFL-numbers and time-steps listed in Table [8] have been measured in this 
regime. For the first test scenario this corresponds only to the first 80 scrt 



(note the slope beginning after about 100 scrt in Figure 10). 



In Table [8] we compare the largest possible time steps Atmax, the average 
time steps Atmean, the maximal and average CFL numbers resulting from the 
adaptive step selection, and the initial CFL numbers. The tests were performed 
for Prandtl numbers Pr = 0.1 and Pr = 0.5 distinguishing Simulations 1 and 
2, respectively, a Lewis number Le = 0.1, Rp = 1.1, and a modified Rayleigh 
number Ra* = 160000 related to the Rayleigh number through Ra* = Ra ■ Pr. 



Comparing the performance of the second order schemes it is obvious that the 
modification of 7 in IMEX SSP2 (2,2,2) has a stunning effect on the stability 
of the scheme. The positive definiteness of dissipation in this method proves 
most effective in suppressing oscillations, permitting an average time-step and 
CFL-number up to two thirds higher than the original IMEX SSP2(2,2,2) 
method. 

Table [s] shows that IMEX SSP2(3,3,2) permits a time-step and CFL number 
more than twice as high as IMEX SSP2(2,2,2) even with modified 7. However, 
a comparison of the computation time given in Table [9] shows that this does 
not improve by the same factor as the CFL-number, since as the time-step 
At grows, the iterative solver for the generalized elliptic problems requires 
more iterations to converge, resulting in an increase in computation time. A 
comparison of the computation times shows that, although IMEX SSP2(3,3,2) 
permits impressively large time steps, the need to solve three additional gen- 
eralized Poisson problems related to the third stage takes its toll, whence the 
method's performance is inferior to IMEX SSP2(2,2,2) with modified 7 and 
in case of Simulation 2, is not even competitive to SSPRK(3,2), though it still 
performs better than the SSPRK(2,2) scheme. Note that for Pr = 0.1, the 
IMEX methods outperform even the best explicit method, while for the more 
moderate Pr = 0.5, the best explicit integrator SSPRK(3,2) is slightly more 
efficient whereas the classical methods noticeably lag behind. 

Interestingly, the initial (preset) CFL number has a negligible influence on the 
actual CFL number reached in the diffusive part of the simulation. However, 
as soon as the fluid velocity seriously restricts At, a higher initial CFL num- 
ber leads to a significantly larger average time-step and reduces the required 
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Method 


CFLstart 


computation time 


number of time steps 


Singlclaycr Pr = 0.1, Le = 0.1, Rp = 1.1, Ra* = 160000 


SSPRK(2,2) 


0.2 


7:17:32 


287 788 


SSPRK(3,2) 


0.5 


4:40:24 


115 691 


IMEX SSP2(2,2,2) 


0.2 


6:04:54 


92 492 


IMEX SSP2(2,2,2), 7=0.24 


0.2 


4:41:15 


63 586 


IMEX SSP2(2,2,2), 7=0.24 


0.3 


4:18:59 


56 741 


IMEX SSP2(2,2,2), 7=0.24 


0.4 


4:10:26 


54 015 


IMEX SSP2(2,2,2), 7=0.24 


0.5 


4:12:38 


53 941 


IMEX SSP2(3,3,2) 


0.4 


4:31:12 


27 673 


IMEX SSP2(3,3,2) 


0.5 


4:27:46 


24 266 


SSPRK(3,3) 


0.2 


10:55:40 


288 607 


IMEX SSP3(3,3,3) 


0.2 


8:19:39 


104 789 


Singlelayer Pr = 0.5, Le = 0.1, Rp = 1.1, Ra* = 160000 


SSPRK(2,2) 


0.2 


7:01:44 


286 Oil 


SSPRK(3,2) 


0.5 


4:34:35 


114 409 


IMEX SSP2(2,2,2) 


0.2 


5:10:01 


75 384 


IMEX SSP2(2,2,2), 7=0.24 


0.2 


4:45:36 


66 997 


IMEX SSP2(2,2,2), 7=0.24 


0.3 


4:42:00 


68 591 


IMEX SSP2(2,2,2), 7=0.24 


0.4 


4:42:24 


66 847 


IMEX SSP2(2,2,2), 7=0.24 


0.5 


4:50:52 


66 828 


IMEX SSP2(3,3,2) 


0.4 


4:43:27 


31 225 


BIEX SSP2(3.3.2) 


0.5 


1: 13:20 


31 112 


SSPRK(3,3) 


0.2 


10:45:43 


286 006 


IMEX SSP3(3,3,3) 


0.2 


7:27:01 


108 852 



Table 9 



Comparisons: computation times and overall number of time steps over 200 scrt. 

computation time, since its value is used to define the time-step restriction for 
the terms integrated with the exphcit part of the IMEX scheme. 

Once the time-step is Umited by t^u\Ai the remaining gain in At by the IMEX 
schemes in ANTARES is essentially due to the optimization of the time-steps 
by the algorithm explained further above. Since At is rather small in that case 
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SSPRK (2,2) 
IMEX SSP2(2,2,2) 
IMEX SSP2{2,2,2), 7=0.24 




100 150 
serf 

(a) CFLstart = 0.2 




SSPRK(3,2) 
IMEX SSP2{2,2,2), Y=0.24 
IMEX SSP2(3,3,2) 



(b) CFL,tart = 0.5 




100 150 
scrt 

(C) CFLstart = 0.2 



Fig. 10. Time-step evolution over 200 scrt in Simulation 1 (see text for definitions). 
Pictures (a) and (b) compare the time-step At of the second order schemes whereas 
(c) shows the evolution of At using SSPRK(3,3) and IMEX SSP3(3,3,3). 
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SSPRK (2,2) 
IMEX SSP2(2,2,2) 
. . , - . , . IMEX gS.P2{2,2,2), 7=0.24 




100 150 
scrt 

(a) CFLstart = 0.2 




SSPRK(3,2) 
IMEX SSP2{2,2,2), Y=0.24 
IMEX SSP2(3,3,2) 



(b) CFL,tart = 0.5 




100 150 
scrt 

(C) CFLstart = 0.2 



Fig. 11. Time-step evolution over 200 scrt in Simulation 2 (see text for definitions). 
Pictures (a) and (b) compare the time-step At of the second order schemes whereas 
(c) shows the evolution of At using SSPRK(3,3) and IMEX SSP3(3,3,3). 
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the convergence of the generahzed Poisson solver is fast enough to allow the 
IMEX schemes to lag only slightly behind their explicit counterparts. 



We point out that optimization of the solver for the generalized Poisson equa- 
tion (21) has the potential for a further significant decrease of computation 
time required, both in the diffusive regime, but also if the time-step is limited 
by Tfjuid- This is important since changing the time integration method during 
a simulation run is not advisable if both the diffusive and the convective phase 
should be interpreted in a consistent manner. 



A comparison of SSPRK(3,3) and IMEX SSP3(3,3,3) also shows that the 
larger time-steps of the semi-implicit method lead to a significant gain in 
computational efficiency in case of a third order method. 



We have investigated the accuracy of the time integration with IMEX schemes 
by a comparison to a reference solution obtained with the SSPRK(2,2) 
method. The reference solution was computed on the same spatial grid but 
with a time-step eight times smaller than that one mentioned in Tables [8] 
and |9] for this method, i.e. for a CFL-number of 0.025. Figure 12 shows the 
root mean square difference of the mass ratio He / (H + He), obtained by 
summation over all grid points and normalization relative to their number, 
between the numerical solutions of the second order SSP IMEX methods and 
the reference solution for each case. For the IMEX SSP2(2,2,2) method both 
the results for the standard choice of 7 = 1 — l/\/2 and the best performing 
value of 0.24 are displayed. We also show the normalized root mean square 
differences between the reference solution and the SSP second order explicit 
methods computed with their standard CFL number given in Table [8] 



For both Simulation 1 and 2 one can easily spot the initial increase of the 
error due to the growing time step for the IMEX methods induced by the 
automatic time step control. A plateau is reached once the time-step stabilizes 
around a typical mean value (cf. also Figures 10 and 11). Note that the IMEX 
SSP2(2,2,2) method with 7 = 0.24 has a smaller error than the original IMEX 
SSP2(2,2,2) method, as expected from the error constant shown in Figure |4j 
The largest differences occur for the IMEX SSP2(3,3,2) method which also has 
the largest mean time-step. The error constants of the different methods (see 
also the summarizing Table 10 below) provide a rough measure for comparing 
simulation runs with similar time-steps. 



However, a comparison for a given point in time has only limited meaning. 
One of the reasons is that the solution changes its nature as a function of 
time. Initial vertical oscillations are damped out (at least first few scrt), then 
the velocity field slowly starts building up (visible after 15 scrt), followed 
by the formation of large scale gravity waves (oscillatory behaviour of the 
error in the range between 25 and 100 scrt) until the waves start to break 
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(a) Simulation 1, first 40 scrt 



(b) Simulation 1, entire run 
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IMEX SSP2(2,2,2) ■ 
... IMEX SSP2(2,2,2) 7=0.24 

'ZX IMEX SSP2(3,3,2} ■ 




5 10 15 20 25 30 35 40 




SSPRK(2,2) 
SSPRK(3,2) 
IMEX SSP2(2,2,2) 
IMEX SSP2{2,2,2) y=0.24 
IMEX SSP2(3,3,2) 



(c) Simulation 2, first 40 scrt 



(d) Simulation 2, entire run 



Fig. 12. Time development of the root mean square difference per grid point of the 
mass ratio He / (H + He) between a reference solution with the SSPRK(2,2) method 
and very small time-step (CFL-number 0.025) and various explicit and IMEX SSP 
methods of second order. In the top row, picture (a) displays the first 40 scrt on a 
linear scale for Simulation 1 and picture (b) shows the results for the entire run on 
a logarithmic scale. In the bottom row, pictures (c) and (d) show equivalent results 
for the case of Simulation 2. 



and turbulence sets in. The importance of the contributions of each of the 
dynamical equations also changes during this development. The whole process 
leads to an increasing error until a statistically stationary, turbulent state is 
reached. For Simulation 1 this occurs at around 150 scrt. For Simulation 2 this 
is just about to occur shortly after the end of the simulation time of 200 scrt 
(the delay of the time development in this case is caused by the larger viscosity 
of the fluid which follows from the choice of Pr and Ra*). In the turbulent state 
of the system spatial correlation is lost on very short timescales (a few scrt). 
Thus, also the reference solution no longer has a meaning due to the chaotic 
behaviour of the solution. The spread of the error and also the saturation value 



observed during this phase (picture (b) of Figure 12) is set by the Dirichlet 
vertical boundary conditions on the concentration c. The third order methods 
IMEX SSP3(3,3,3) and SSPRK(3,3) behave analogously. 
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(a) SSPRK(2,2) 



mfhe 
0.7303 




(b) IMEX SSP2(2,2,2) 




(c) SSPRK(3,2) 



(d) IMEX SSP2(3,3,2) 




(e) SSPRK(3,3) (f) IMEX SSP3(3,3,3) 

Fig. 13. Simulation 1 at i = 125 scrt. 

To further illustrate that the large time steps of the IMEX methods during 
the diffusive phase do not degrade the accuracy of the time development of 
the solution, we compare the simulation results obtained by the different time 



integrators in Figures 13 and 14 The pictures show the mass ratio of He vs. 



He+H at each spatial point at a given instant in time. Figure 13 demonstrates 
that after the end of the diffusive phase, just at the onset of turbulence, which 
occurs at ~ 125 scrt for this problem, the results are quite comparable. The 
most visible differences can be found for the case of the simulation with the 
largest time-steps (picture (d)). It also shows that the root mean square errors 



displayed in Figure 12 are negligible on a qualitative (and rough quantitative) 
level as long as they are smaller than about 10"'^. As expected, however, some 
time after the onset of turbulence the solutions necessarily have drifted apart. 



This is demonstrated in Figure 14, where the solutions already look different 
from each other. 



Recalling Figure 12 it is not surprising to find the largest differences in Fig- 
ure 13 for the schemes having the largest time-steps during the diffusive phase. 



Still, the large scale structures of the solution begin to diverge only once the 
turbulent phase has been reached which in turn for each of the different time 



31 



1^ 




(a) SSPRK(2,2) 



(b) IMEX SSP2(2,2,2) 





mfhe 
0.7303 



(c) SSPRK(3,2) 





(d) IMEX SSP2(3,3,2) 




(e) SSPRK(3,3) 



(f) IMEX SSP3(3,3,3) 



Fig. 14. Simulation results corresponding to those in Figure 13 at t = 200 scrt 



integration methods occurs after about the same integration time t. Large 
time-steps during the diffusive phase hence yield acceptable accuracy. Indeed, 
the spatial resolution is more important than the temporal one. This can be 
demonstrated by using a high resolution grid of 800x800 points. We have per- 
formed such a reference run for the case of Simulation 1 with the SSPRK(2,2) 
method for time integration. Note that the doubling of resolution leads to a 
four times smaller time-step during the diffusive phase. Looking at the first 



row (pictures (a) and (b) of Figure 15) the differences at 100 scrt, i.e. just 
after the onset of wave-breaking and the beginning of the turbulent phase, are 
still small: in the simulation with higher spatial resolution the breaking tips 
are somewhat more pronounced and the contrasts sharper. This has changed 



at 125 scrt shown in the second row of Figure [T5| The common initial condi- 
tion may still be inferred, but the simulations have notably evolved away from 
each other. Clearly, the spatial resolution is much more important than the 
influence of the time steps and the time integration method chosen, since pic- 
ture (c) of Figure 15 is essentially indistinguishable from its counterpart with 



standard time resolution, picture (a) of Figure 13 Furthermore, at 100 scrt 
the simulation using the IMEX2(3,3,2) method for time integration, which 
has the largest time-step during the diffusive phase, is nearly indistinguish- 
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(a) high time resolution, t = 100 scrt 



(b) high spatial resolution, t = 100 scrt 



mfhe 
,7303 




0,73t 
I0,7 




(c) high time resolution, t = 125 scrt (d) high spatial resolution, t = 125 scrt 




i,730: 
0,7 




(e) high time resolution, t = 200 scrt (f) high spatial resolution, t = 200 scrt 

Fig. 15. Simulation 1 with SSPRK(2,2) time integration and high temporal resolu- 
tion (left column) as well as high spatial (and temporal) resolution (right column). 
The three different rows show the results at different time t in units of scrt. 

able from the SSPRK(2,2) run with high time resolution shown in picture (a) 



of Figure 15 (the IMEX results are not shown here for the case of 100 scrt, 



since the differences to the latter picture are very difficult to spot). 



We conclude that the spatial resolution is indeed more important than high 
temporal resolution and large time steps during the diffusive phase are clearly 
tolerable for simulations of astrophysical convective flows, if they do not af- 
fect stability. Resolutions of 800 grid points per spatial direction in 3D are 
usually not affordable anyway and sometimes even large parametric studies 
in 2D may still be too expensive (for instance, for the case of semi-convection 
and purely explicit time integration methods). We note that the necessity of a 
resolution of 400 points, which has been used for most of the simulation runs 
shown here, was calculated following |38] and [19] where in turn the physical 
arguments of [43j had been used to estimate the thickness of solute and ther- 
mal boundary layers in semi-convection and the applicability of this approach 
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to the parameter range we are interested in had been confirmed. Thus, for the 
parameters of the more demanding one of our models, Simulation 1, we con- 
cluded the smallest structures of interest, the solutal boundary layers, to span 
6 grid points, if the whole box is discretized by 400 points in each direction. 
Examples for them are the top and bottom boundary layers which can easily 



be seen in Figures [14] and [15] Indeed, in the simulation with a high spatial 
resolution of 800 points in each direction these layers are hardly any thinner 
than in the case of 400 points (cf. the bottom row of Figure [Ts] which shows 



the simulations developed well into their final, turbulent state). 

The duration of the diffusive phase is determined by the stability of the strat- 
ification, parametrized through Rp, and related to the buoyancy term, the 



second term on the right-hand side of (15) and thus also the second term of 
F{y{t)) in the same equation. Timescales related to this term can be computed 
for a variety of physical problems. They include the reciprocal of the growth 
rate of small density perturbations when a fluid of higher density p2 is layered 
above fluid of lower density pi (this situation is denoted as Rayleigh- Taylor 
instability, see [6]). In this case the growth rate follows the dispersion relation 
= gk{p2 — Pi)/(p2 + Pi) for a local gravitational acceleration g and a per- 
turbation with wave number k. Its magnitude can be bounded by the simple 
relation u"^ = gk. Similar dispersion relations are found for gravity waves and 
growth rates of convective instability (see the classical paper [7J and also [26] 
for a summary). As implied by u'^ = gk and the form of (15), one can esti- 



mate a buoyancy time-step restriction by tbuoy = min{(Aa;)^}/(7^/^ (see also 
[52]). For the present simulations we find tbuoy ~ 360 s, i.e. about 0.069 scrt. 
Note that this is about four times larger than the largest time-step reported 
in Table |8] Though irrelevant in the asymptotic limit and not a constraining 
quantity here, one might still consider this term to be integrated implicitly 
in other cases. However, if the excitation and breaking of waves is important 
to describe the onset of the turbulent convective flow, as in the present case, 
damping of such waves by an implicit time integration may be undesired. 
Hence, explicit time integration of the buoyancy term could be preferred for 
physical reasons, even if the time-step were actually constrained by such a 
splitting for the time integration. 



Numerical Results with non-SSP Schemes 



From Tables [8] and [9] one can readily see that during the diffusion domi- 
nated phase the SSP IMEX schemes achieve time-steps which exceed the re- 
gion of absolute monotonicity ensured by (32) for IMEX SSP2(2,2,2) and 



their counterparts given after ([33]) and ([35]) for IMEX SSP2 (3,3,2) and IMEX 
SSP3(3,3,3), respectively. One might hence question whether the property of 
absolute monotonicity is really necessary for the time integration of the nu- 
merical simulations we have considered above. To show that this property is 
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indeed required we have performed several test runs for the case of Simula- 
tion 1 with time integration schemes which do not diminish the total-variation 
norm. 



The first candidate we have investigated is the ARK3(2)4L[2]SA scheme pro- 
posed in [25]. This is an L-stable, stiffly accurate third order, additive Runge- 
Kutta method with four stages. With its choice of coefficients it belongs to the 
group of IMEX methods. However, neither its explicit nor its implicit part are 
strong-stability-preserving (the Butcher arrays feature negative coefficients, 
cf. Theorem 4.2 in EHI), hence also the entire method does not fulfill the crite- 



ria of Theorem lA on absolute monotonicity. If the SSP-property were of no 
importance, this method should be quite robust. However, it turns out that 
this is not the case when we use it to integrate Simulation 1 in time. Tak- 
ing CFLstart to be 0.2 as for the other IMEX methods presented in Table |9} 
the time integration with ARK3(2)4L[2]SA crashes after just 78 time-steps. 
Indeed, CFLstart has to be lowered to 0.1 to successfully launch the simula- 
tion. But over the first 10 scrt CFLmax is found to never exceed ~ 0.2 and 
the average CFLmean is only ~ 0.15. In conclusion the size of the time-steps 
achieved with this method have been found to not exceed those achieved with 
the explicit, three-stage, third-order SSPRK(3,3) method of [S]. Compared 
to IMEX SSP2(3,3,2), the maximum and mean CFL numbers are 25 and 20 
times smaller, respectively. We have thus given up this simulation run after 
10 scrt: evidently, the ARK3(2)4L[2]SA method is not efficient for the kind 
of problems we are interested in. The implicit, L-stable and stiffly accurate 
nature of this scheme is not sufficient to provide any advantages on its own 
during the diffusive phase of the simulation. 



To further investigate the importance of the strong-stability-preserving prop- 
erty we selected the classical, explicit, third order, three-stage Runge-Kutta 
method first proposed in [18] and known as Heun's third order method. This 
is a non-SSP scheme since not all of its stages are used in the final integration 
which yields ?/new, as pointed out in [28], where it was used to illustrate the 
growth of solutions measured in standard norms for both parabolic and hyper- 
bolic problems in situations where the exact solution is not growing in these 
norms. By comparison the SSPRK(3,3) scheme was found to not exhibit such 
growth. When we apply Heun's third order scheme to integrate Simulation 1, 
we achieve a stable simulation over the entire extent of 200 scrt with the same 
average time-step and with the same Courant number as for the SSPRK(3,3) 
scheme. However, during the diffusive phase in Simulation 1 the solution it- 
self is slowly growing in time while a velocity field is being built up by the 
convective instability and at the grid scale the dissipation is provided by the 
parabolic terms during the entire simulation. 

The semi-convection problem discussed here is a rather benign example for nu- 
merical simulations in astrophysical applications. If we apply the same scheme 
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to a simulation of solar surface convection as in [35], i.e. for a case of 219 x 159 
grid points and a standard, moderately low resolution of 18.57 x 40 km^, and 
a standard choice for the microphysics with non-grey radiative transfer for 
the calculation of Qrad, the differences in stability become apparent. For this 
physical problem the term representing viscous dissipation in the momentum 
equation does not provide sufficient dissipation on the grid-scale at any resolu- 
tion achievable in the foreseeable future and the dissipation properties of the 
temporal and spatial discretization of the advection operator become impor- 
tant (the term implicit large eddy simulation is used in such cases). While the 
SSPRK(3,3) method and also the SSPRK(2,2) method, used together with 
the spatial discretization of [IQ], have no problems in completing a simulation 
of 20 scrt with a CFL number of 0.25, Heun's third order method leads to a 
crash after 9.2 scrtj^Such failures are usually caused by numerically induced 
fluctuations in the low density and temperature region of the simulation box 
which result in negative or at least unphysically small values, whence they 
fall outside the tabulated region of microphysical properties. We consider this 
flnding as sufficient to exclude non-SSP methods from being recommendable 
for numerical simulations of stellar convection, since also here we have chosen a 
rather benign test case within its class. Numerical simulations of stellar surface 
convection in white dwarfs, A-type stars, or Cepheids reach far more extreme 
conditions with up to three times higher, super-sonic Mach numbers, density 
contrasts around shock fronts higher by an order of magnitude and more, and 
for the case of A-stars and Cepheids, at lower effective resolution because of 
four to ten times steeper gradients and limited computational resources. 



We finish our study of non-SSP Runge-Kutta schemes with a third case: an 
IMEX Runge-Kutta method where both the explicit and the implicit part 
are strong-stability-preserving, but the combined scheme does not fulfill the 
criteria for absolute monotonicity in the sense of Theorem |1.1[ The scheme 
was proposed in [M] (Table IV, page 139) and indeed the IMEX SSP2(3,3,2) 
scheme (33) is a modification of that scheme proposed in [21] to obtain a 
nontrivial region of absolute monotonicity. The original scheme of [36] differs 
from that one only by having the entries {1/4, 0, 0} and {0, 1/4, 0} in the first 
two rows of the Butcher array of the implicit scheme instead of {1/5, 0, 0} and 
{1/10,1/5,0}, respectively. This scheme performs quite well. Indeed, during 
the diffusion dominated phase of Simulation 1 its time-steps are even 6% to 
8% larger than those achieved with IMEX SSP2(3,3,2), while during the tur- 
bulent phase they fall back to at most the size achieved by the scheme (33). 
Once more, however, we recall that astrophysical simulations often have to 
deal with limited resolution at least in part of the simulation domain. To re- 
produce such a case we have run Simulation 1 with both (33) and the original 



^ We would like to thank H. Grimm-Strele for performing the 2D solar convec- 
tion simulations with ANTARES to test the stability properties of Heun's non-SSP 
explicit third order Runge-Kutta scheme. 
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scheme of [36] for the case of only 100 x 100 grid points while leaving every- 
thing else unchanged. In this case the boundary layer due to concentration 
c is represented vertically only by one to two grid points (see the discussion 
on resolution and reference solutions further above). While during the diffu- 
sive phase no major differences become apparent, the behaviour found for the 
advection-dominated, turbulent phase was discrepant: whereas nothing suspi- 



cious occurred for the scheme (33), the time-step of the scheme of |36] dropped 
to arbitrarily small values after ~ 113 scrt, which indicates the occurrence of 
two-point instabilities, and the simulation had to be terminated. 

We conclude that only Runge-Kutta methods which are strong-stability- 
preserving have the necessary prerequisites for stable time integrations of as- 
trophysical convection simulations. If, in addition, the time-step is limited 
by diffusion processes, this limitation can be overcome by IMEX methods 
provided their explicit and implicit parts are strong-stability-preserving. To 
ensure stability also in cases of low resolution IMEX SSP methods should also 



have a nontrivial region of absolute monotonicity as required by Theorem 1.1 



5 Conclusions and Outlook 



In this paper we have given an extensive discussion of the mathematical prop- 
erties and practical usefulness of total-variation-diminishing implicit-explicit 
Runge-Kutta methods for the time integration of advection-diffusion equa- 
tions arising in the simulation of double-diffusive convection in astrophysics. 
In this section, we summarize the results obtained in Sections [3] and |4] (stabil- 
ity, dissipativity, accuracy and efficiency), and give a brief outlook on future 
developments. 

The stability regions for the IMEX methods, their implicit sub-parts and the 



explicit schemes are given for comparison in Figure 16 The left boundaries 
of the stability regions zi^ft, the points where the amplification factors from 
the dissipativity analysis become zero and their moduli exceed 1, and the 
error constants C for all the methods we have investigated are summed up in 
Table M 

We found that methods introduced in [20|2H22||36] excel over the classical ex- 
plicit methods [4iJ. It was found that among explicit schemes, only the explicit 
SSPRK(3,2) scheme first proposed in [2S] is competitive in situations where 
explicit time integration can be expected to yield sufficient efficiency and ac- 
curacy. Examples for such a scenario include simulations of solar granulation 
at moderately high spatial resolution, where the time-step limitation associ- 
ated with diffusion is negligible (see [S^ for results on this problem obtained 
with the ANTARES code and [H] for a review on the underlying physics). 
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(d) A in (29) 



(g) SSPRK(2,2) 




(a) SSP2(2,2,2) (29) (b) SSP3(3, 3, 2) (33) (c) SSP3(3, 3, 3) (35 ) 



(e) A in (33) 



(f) A in (35) 



(h) SSPRK(3,2) 



(i) SSPRK(3,3) 



Fig. 16. Stability regions for IMEX methods (first row) and their impHcit and ex- 
pHcit sub-parts (second and third row). 



From Figure 16 and Table 10, clearly there is no single scheme which features 
the most advantageous properties in all considered aspects. However, we found 
in numerical experiments that the most efficient method seems to be (29) with 
the choice 7 = 0.24. This value deviates from the optimal value for strong 
stability, but leads to a scheme with favourable dissipativity, stability, and 
accuracy properties. Depending on the domain of stability required for a given 
problem the value of 7 in (29 ) may be optimized such that it is sufficiently large 
for stability, but small enough to minimize the error constant, while showing 
favourable dissipation properties (a strictly positive amplification factor with 
modulus less than 1 for any wave number k other than zero). 



We note that for numerical problems arising from a method of lines approach 
to the equations of hydrodynamics, as discussed in this paper, lower order 
methods usually have sufficient efficiency to be competitive, since the spatial 
discretization limits the overall accuracy. Hence, the best explicit scheme we 
have tested for this kind of apphcation is SSPRK(3,2), as it permits the largest 
CFL numbers among methods of this class at an affordable computational 
cost and with sufficient accuracy. By comparison, the classical methods of 
second and third order [H] offer the convenience of being usable together as an 
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embedding formula. However, this approach is more than twice as expensive 
as SSPRK(3,2), as can be seen from a comparison with SSPRK(3,3) using 
Table E 



Method 




54th = 


fl'4th = 1 


C 


IMEX SSP2(2,2,2) 


—00 


0.452* 





5.17 


IMEX SSP2(2,2,2), 7 = 0.24 


-50 


— 


9.375 


2.79 


IMEX SSP2(3,3,2) 


— CXD 


0.455* 




8.05 


IMEX SSP3(3,3,3) 


-3.248 


0.348* 


0.609 


11.6 


Forward Euler 


-2 


0.187* 


0.375 


12.6 


SSPRK(2,2) 


-2 




0.375 


16.2 


SSPRK(3,2) 


-4.519 


0.672* 


0.847 


6.40 


SSPRK(3,3) 


-2.512 


0.299* 


0.471 


22.8 



Table 10 

Summary of the analysis of SSP integrators. The asterisk in the third column in- 
dicates a change of sign at /u for g4th(±vr, /x). Other details are given in the text. 



We have also demonstrated that the larger time-steps achieved by SSP IMEX 
methods reduce the accuracy of the solution during the diffusive phase of the 
semi-convection simulations by an acceptably small amount. For the applica- 
tions shown here, and indeed for a majority of astrophysical fluid dynamical 
simulations, the accuracy is limited by spatial resolution (and thus eventually 
by existing computational resources) while the time-steps are limited by sta- 
bility. This makes IMEX methods attractive, since quite often the most severe 
limitations stem from stiff terms representing diffusion processes (for restric- 
tions due to sound waves other operator splitting based methods are existing). 
However, as we have shown by a comparison with results from non-SSP meth- 
ods, it is important that the IMEX methods are strong-stability-preserving 
to maximize stable time-steps no matter whether the constraint is due to 
the implicitly integrated terms (diffusion) or the explicitly integrated ones 
(advection). From that point of view SSP IMEX methods with a non-trivial 
region of absolute monotonicity as defined by Theorem IT are the most robust 
methods, because they allow achieving optimally large time-steps also at low 
resolution. We note here that while the region of absolute monotonicity has to 
be observed with respect to the explicitly integrated advection operator of the 
dynamical equations, stable time integration can be performed with step sizes 
falling outside it, if the restriction is due to the implicitly integrated diffusion 
terms. Thus, the class of optimal integrators for this kind of problem is prob- 
ably larger than that of the SSP IMEX methods. However, none of the other 
time integration methods could significantly outperform them with respect to 
the time-steps achievable, and we have always found at least one case, where 
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competing methods fell substantially short or even failed. 



There is potential to further optimize the implementation of IMEX methods. 
The additional computational effort due to the implicit subpart is compen- 
sated for by accuracy and stability, but could be reduced in the future by 
replacing the solver for the linear equations associated with the arising gener- 
alized Poisson problem by a multigrid solver. In the Boussinesq approximation, 
additional solution of a Helmholtz equation is necessary instead. This widens 
the choice of fast solvers for the system of linear equations introduced through 
implicit time integration. The benefits expected from faster solvers would al- 
low taking full advantage of the potential of method (33 ) that is implied by the 
large time steps reported in Table [8j Such an improvement would likewise be 
useful for the present problem to minimize the overhead by any of the implicit 
schemes in the regime where the time-step is limited by Tfluid rather than tt- 
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Appendix 



The Explicit SSP Scheme of IMEX SSP2(3,3,2) 

We also give the corresponding results for the explicit SSPRK(3,2) scheme A 



from (33), since in Section |4| we show it to excel in its practical value over 
the classical explicit SSP Runge-Kutta schemes schemes [H]. This scheme 
was first published in [28] and later declared the optimal second order scheme 



with three stages in ^2] as well as in [37], and independently also in [12] . 
The stability function is 



i?^(2;) = l + 2; + -+^2- 
The stability region where |-R(;i;)| < 1 occupies a bounded region in the nega- 



tive half-plane, and is tangent to the imaginary axis, see Figure 17, Note that 
lim3}(^)^_oo \R{z) \ = oo. 
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Fig. 17. Stability region of method SSPRK(3,2) (scheme A from (33)). 



Table 11 



9 


9(0, f^) 





1 


7r 
4 


1 + fiV2-2fi + 3n^ -2 fi^V2 + 7/6 /^^^ - 5/3 


TT 

2 


l-2^ + 2;u2-2/3/i3 


TT 


l-4/x + 8Ai^- 16/3^3 



Values of g{0,n) for some 0, SSPRK(3,2) scheme (explicit scheme A from (33)), 
three point space discretization ([s]). 



e 


ff(^,/^) 





1 


TT 
4 


1 - 5/2 ^ + 4/3 ^/2 + if /i2 - 10/3 V2 - ^ ^3 + If 


TT 

2 


l-7/3^ + i/.2-iiM^ 


TT 


l_16/3^+lf ^2 _ 1024 ^3 



Table 12 

Values of g{6, fx) for some 0, SSPRK(3,2) (explicit scheme A from (p3|)), fourth order 
space discretization ([9]). 

The point where the stabihty region intersects the negative real half-hne is 
located at x ~ —4.519. 



The dissipativity analysis for A yields the amplification factors for the stan- 
dard three-point space discretization ([s]) and the fourth order stencil These 
amplification factors are evaluated at the points 6 G {0, |, |,7r} in Tables 



and 12, respectively. 
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For the three-point space discretization ([8]), the first positive zero of g{'n',n) 
is ~ 0.8968, where the function changes its sign, and 5^(71, /i) = —1 for /i ^ 
1.129. The first positive zero of g{7r,n) is ~ 0.6726 for the fourth-order space 
discretization ([9]), where the function changes its sign, and g{rc,iJ,) = —1 at 
/i ^ 0.8474. 
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